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SUMMARY 

The methods of Hamilton and Lagrange in the calculus of 
variations are employed to optimize planetocentric orbital transfers 
using fixed-thrust engines, where the thrust acceleration << local 
gravitational acceleration. The equations were derived for minimum 
time transfers. However, since fixed-thrust engines are used, the 
minimum fuel problem is solved simultaneously. Numerical solutions are 
obtained by integrating the optimized differential equations of motion 
for the problem, the actual integrations being performed by an IBM 360/30 
computer. The numerical integrations result in optimal spiral 
trajectories representing orbital transfers between specified coplanar, 
co-axial ellipses, with the optimal thrust angle being calculated at all 
points along the trajectories. For an orbital transfer at the planet 
Mars, the low fixed-thrust system yields a significant increase in the 
amount of payload that can be delivered into the desired final orbit 
compared with a transfer accomplished by an impulsive (high thrust) chemical 
rocket propulsion system. However, the time required to achieve the 
transfer is necessarily longer in the low-thrust case. It is then necessary 
to examine the trade-off between the increase in payload and the extra 
transfer time to determine whether or not the fixed-thrust system is 



desirable 
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I. INTRODUCTION 

As we explore Che space around us, we look for ways of doing things 
more efficiently, ways to save money, and ways to put heavier payloads into 
orbit at the same or less cost than is presently possible. One feasible way 
of accomplishing this is through the use of solar-electro-static or nuclear- 
electro-static propulsion systems. They require much less fuel (weight- 
wise) than chemical rocket engines, have much higher specific impulses, and 
can result in considerable increases in final payload delivered. They do, 
however, require more time to achieve the same orbital maneuvers because of 
the much lower thrust inherent in these type of systems. 

Considerable research on optimal low thrust maneuvers has been done 
in recent years. A large part of this work has been done by T. N. Edelbaum. 

He has published papers which contain solutions for optimum power-limited 
transfer and rendezvous between arbitrary coplanar and coaxial elliptic and 
circular orbits. The assumptions involved in the development are 1) the 
thrust acceleration is much smaller than the local gravitational acceleration. 
2) the propulsion system is power - 1 imi ted , i.e., it operates at constant 
exhaust power with a variable thrust magnitude inversely proportional to the 
exhaust velocity. 3) the direction and magnitude of the thrust, both of 
which are assumed to be completely variable, is to be determined as a function 
of time so as to minimize the fuel consumption for a fixed total transfer time. 
Edelbaum has also developed an analytic solution for fixed thrust systems, but 
this is only applicable to transfers between almost circular orbits, for small 
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changes in a and e , and for large changes in a and i (orbital plane 
inclination) (Reference 1, 2, 3). 

This thesis addresses the problem of optimizing fixed, low- 
thrust orbital transfers with large changes in a and e . Based upon a 
Lagrange-Hamilton formulation developed by Dr. M. Handelsman (Reference 4) 
the numerical solution of optimized transfer trajectories is performed by 
computer. The solutions involve the optimization of plane tocentric transfers 
between coplanar, coaxial ellipses or circles of specified characteristics, 
with no time constraint. Because of the relatively small fuel expenditure 
in these transfers, it is possible to replace "fixed thrust 1 ' with "fixed 
thrust acceleration" throughout the entire transfer with little error. Solar 
electric engines lose a certain amount of thrust as their distance from the 
sun increases; this could hardly be considered constant thrust in heliocentric 
space. But, since the problem has been limited to a plane tocentric study, 
the thrust of the solar electric engine can be assumed to be constant within 
the planetocentr ic space. 

By fixed thrust, it is meant that the thrust magnitude is constant 
and that the engine remains turned on at all times during the transfer. The 
problem then remains, where to point the thrust at all times. Now, there are 
many factors which could be optimized in an orbital transfer of this type. 

In this thesis two separate phenomena were optimized. First, orbital trans- 
fers were optimized for minimum time using fixed thrust engines. Since the 
thrust is fixed, and, therefore, continuous, then a minimum time transfer 
would also be a minimum fuel transfer, the amount of fuel expended being 
directly proportional to the length of time the engine is turned on. Secondly, 
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orbital transfers were optimized for minimum E (eccentric anomaly); in 
other words, transfers were achieved to the desired orbit in the least amount 
of revolutions . The significance of this type of transfer trajectory is that 
it lends considerable insight into the minimum time, minimum fuel problem. 

The problem formulation involves deriving the perturbation equations 
of motion for fixed thrust engines of low, constant magnitude, and then opti- 
mizing them for minimum time and minimum E using Hamilton > Lagrange, and 
Mayer techniques in the calculus of variations. (Reference 4). After being 
optimized for minimum time and minimum E respectively, the differential 
equations are set up to solve for the optimum thrust angle necessary at all 
times to achieve the optimized transfers. These sets of equations are then solved 
by numerical integration, in this case on an IBM 360-30 computer system. An 
analytic solution is not available at present because of the inability to solve 
for Langrangian multipliers. This is the reason for the numerical approach. 

The basic assumptions follows; 

1) The thrust acceleration is much less than the local gravitational acceler- 
ation. (A«g). 

2) Continuous, fixed-magnitude thrust. 

3) Relatively small total fuel consumption; therefore, for the sake of simpli- 
city, the thrust acceleration is assumed to be essentially constant over the 
transfer. That is, constant thrust magnitude is taken as synonymous with 
constant thrust acceleration. (A program modification to incorporate the 
changes in thrust acceleration due to fuel depletion is simple to implement). 

A description of the problem and its formulation follows. 
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II. DESCRIPTION OF PROBLEM 

An initial elliptic orbit is assumed about a planet or any other 
body which has an inverse-square gravitational field. The ellipse has 
specifications (a Q , e^) , where a^ is the initial semi-major axis and 
e^ is the initial eccentricity. The problem then is to transfer optimally 
to a desired final coaxial, coplanar ellipse with specifications (a., e^) . 

See figure 1. This is accomplished by turning on the engine at the initial 
ellipse periapse and letting the optimally-directed fixed thrust program guide 
the vehicle along a spiral trajectory until the final ellipse is obtained. 

The engine is then shut off and the spacecraft is in its new orbit, the trans- 
fer having been conducted in minimum time and fuel expenditure or in the least 
amount of revolutions, depending upon which program is used. The numerical 
results show that the minimum time transfers always take more revolutions 
than the minimum E transfers, and the minimum E transfers always take more 
time than the minimum time transfers. This is as expected and lends consider- 
able support to the calculations. 

Each particular transfer is determined by the selection of the initial 
Lagrangian multiplers. Given an initial ellipse, different sets of inital 
Lagrangian multipliers will cause the vehicle to transfer to different ellipses. 
Since we have no analytic way of picking these multipliers at present, many 
cases are run to cover a wide range of transfers from given initial orbits. 

The formulation of both the minimum time and minimum E equations 
follows as the first phase of the problem solution. 



SPIRAL TRANSFER 
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FIGURE 1 
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III. LAGRANGE FORMULATION IN VARIABLE E 
A. Orbital Perturbation Equations 

The general orbital perturbation equations in terms of radial 
tangential thrust forces, S and T respectively, shown in Figure 2, 
are as follows: (Reference 5) 

Some useful relations follow: 




FIGURE 2. 

Vehicle Force Vectors 



a = 



da 

dt 



(S e Sin f + T -Jr ) 






e = 



de _ 

dt na 



S s in f + T(cos E + cos f) 



Us ing: 



(i) 



( 2 ) 




cos E + 



cos f 



e 



( 



r 
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p = a(l-e ) 



f = true anomaly = 9 



E = eccentric anomaly 



Equations (1) and (2) become: 



2e|sin f| S 2a V l-e^ 

a = — — — 1— - _L ■ 



a = 



n yJi-( 



r n 



(3) 



\C 2 



sin f 



na 



S + 



VT7 2 

2 

a ner 



2 /i 2 n 2 

a (1-e ) - r 



(M 



2 3 

na =y=/u=,k = GM 



(3) 



n = 



V k/a^ 



= 2 7T / 



(period) 



S F NORMAL/ m F sin« / 



( 6 ) 

(7) 



T = F TANG./ m _ F cos« / m 



( 8 ) 



For an ellipse: 



r = a (1-e ) 

1+e cos f 

Through substitution and simplification equation (3) yields: 



(9) 



a 

a 



2, x JX 

Vl-e 2 m k 



e sin 9 + sin« + (1 + e cos 9) cos« 



( 10 ) 
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Similarly, equation (4) yields: 



e - \[T- ? -Z-y/l 

m * k 



^ ,2 cos 0 + e(l + cos 0) cos a 

Sin 0 s in a + ; — * L 

1 + e cos 0 



( 11 ) 



Next, it is desired to write equations (10) and (11) in terms 
of the eccentric anomaly E. The following relations between E and 0 
are used (Reference 6): 



1 + e cos 0 

sm 0 = — ?==-t sin E 

V 1“ e 



( 12 ) 



1 + e cos 0 = 



1-e 



s in 0 = 



cos 0 = 



r = 



1+e cos 0 



1 + e cos 0 



1-e 


cos E 


Vl-e 2 


1 

sin E 


1-e 


cos E 


cos 


E - e 


1-e 


cos E 




Ml-e 2 L 


9 


1-e 2 




1-e cos E 


1 - 


2 

e 


1 - 


e cos E 



= a(l-e cos E) 



(13) 



(14) 



(15) 

(16) 



(17) 



sin E = sin 9 



dQ = 



x/TI 2 



1-e cos E 



dE 



(18) 

(19) 



I 



I 
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B. Comparison with Edelbaum's Equations 

Substituting equations (12-19), equations (10) and (11) can be 
written in terms of E: 



a. 

a 



2F 



(1-e cos E) m 



e sin a sin E + cosor 



,\C 2 



( 20 ) 



e = 



(1-e cos E)m 



2 - • - • . + \/r? - - 2 



(1-e ) sin E sin 



(2 cos E - e cos E-e) cos a 

( 21 )- 



Equations (20) and (21) agree with Edelbaum's equations for 
low thrust perturbation (Reference 1). 

Since 






dt = dE i (1-e cos E) 



( 22 ) 



equations (20) and (21) become; 



da 

dE 



2Aa' 



e sin E sin 



,r 2 

Lna + \l-e cos a 



(23) 



de _ Aa 
dE k 



2 r 



2 / 2 2 

(1-e ) sin E sin a + Vl-e (2 cos E-e cos E-e) cos a 



(24) 



Since we are optimizing the thrust angle for minimum flight time 
the following integral is minimized: 

E 4 



^ Wa 3 / ^ (1-e cos E) dE 

V 



« 



(23) 
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This is of the form (Reference 4, 7) 



t 



f 




X 2 , E) dE 



(26) 



where : 



f 

o 



5 P 




(1-e cos E) 



(27) 



X x = a; X 2 = e (28) 

C. Employing the Hamiltonian 

To find the minimum extremal, we write the Hamiltonian for this system of 
equations , 




X. 



i 



(E) f. 



(29) 



where 



f 



1 




da 

dE 



(30) 



f 



2 




de 

dE 



(31) 



The independent variable is the eccentric anomaly E. Substituting 



into II: 

I 3 da de 

= - \a /k (1-e cos E) + * a dl + X e dE 



H 



(32) 



I 



i 
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0 ^ da , de 

Substituting for — and — 

dL dE 






H = - ^ aJ / k ( 1_e cos E) + X a 



3 r 



e sin E sin a -I- 






1-e COS or 



+ A 



Aa 



e k 



(1-e ) sin E sin a + 



V 



2 2 
1-e (2 cos E - e cos E - e) cos a 



(33) 



D. Euler-Lagrange Equations for Adjoint Variables 



The following are the Euler-Lagrange equations for the adjoint variables \ 



d A 4 

] 

dE 



dH 

ax. 

l 



(i = 1 , 2 ) 



(34) 



Expanding equation (34): 



d A , 
dE 



df -i 



d*. 



( A " 3 — + A 7 

a a a e da 



) + — ( X f. + \ f 9 ) 



d& 



al A e 2 ' da 



(35) 



d X 

€ 

dE 



= - ( X 



d«, 



0 f 



+ X 



2,1 



) +“ ( A f, + A f 0 ) 



dP 



ade A e 3 e f} A a 1 e 2' 3 e 



(36) 



E. Transversality Condition 



Next, the transversality condition. 



HdE + A . dX. 

l l 



= 0 



(37) 
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which is to be evaluated using the boundary conditions: 

t = 0: a = a , e = e , E = 0 (38) 

o o o v ' 

t = t f : a = a , e = e , E f = OPEN (39) 



The choice E q = 0 is made for convenience; this is the periapse point, the 
usual choice for the initial retro- impulse into a bound planetocen tr ic elliptic 
orbit from an incoming hyperbolic trajectory. 

Combining equations (37), (38), and (39), 

A a f l + A e f 2 " P = £ = 0 ( at t=t f) ( 40 ) 

F. Control Equation for Optimal Thrust Angle 

Equation (40) is a necessary condition for a transfer to be optimized for 
minimum time (and fuel). For the optimal thrust angle a , 



3 H 
3 « 



0 



3 f 3^2 

A a +A e 3« 



(41) 



which becomes 
A a 2Afl3 



A ^a 2 I 2 1 2 

[(l-e ) sin E cos a - Jl-e (2 cos E - e cos E - e) sin a] (42) 



[e sin E cos a 



V 



. 2* . 
l-e sin 



t 



k 
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Divide through by cos or and solve for tan a : 



tan a = 



2a A e sin E + A (1 -e ) sin E 

3 , 6 



* a yjl-e ^ ^ e ^1-e (2 cos E - e cos^ E - e) 



( 43 ) 



To simplify the equation for computer programming purposes let: 



X 


— e sin E 


(44) 


z 


_ 2 
— (1-e ) sin E 


(45) 


Y 


CM 

<D 

L i 

ill 


(46) 


W 


1 2 2 
3 ^(l-e (2 cos E - e cos E - e) 


(47) 


N 


r 2a A X + A Z 

a e 


(48) 


D 


E 2a A Y + A W 
a e 


(49) 


Substituting into equation (41): 




2a A X + A Z M 

a e — N 

tan a - 2a A Y + A W D 

a e 


(50) 



This is the control equation for the optimal thrust angle « . The step variable 
of integration is E. Time of transfer (t) is determined by integration of 
equation (22) which can be written; 



d_t- _ _1_ 
dE p 



( 51 ) 









i 

I 



' 
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Thus, we now have sufficient equations for numerical solutions for 

different transfers, depending on the values assigned to \ and \ 

a 

o 

The problem is next solved using the Mayer formulation to develop an 
alternate but equivalent set of equations. 
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IV. MAYER FORMULATION IN VARIABLE t 



The Mayer formulation of the problem proceeds as follows: 



dX. 

ir - W V 



( 52 ) 



Whe re : 



X L = a, X 2 = e, X 3 = E, = a 



( 53 ) 



The independent variable now is t rather than E as in the 
Lagrange formulation; here we minimize the value of t^, which will be 
called G. Thus: 



G = t . 



( 54 ) 



A. Hamiltonian and Euler-Lagrange Equations 

The Hamiltonian is different and does not contain f 



H = £ 



i=l 



X .(t) f . 
l i 



( 55 ) 



The Euler-Lagrange equations for the adjoint variables \ are 



dA. 

i 

dt 



9 H 
9X. 



( 56 ) 



Thus : 

d *a 

dt 



P 



^ ^ a da + ^e 0a 



) + — o ( ^ f i + A f 0 + A ) ^ 

7 ^2 ' a 1 e2 E d a 



( 57 ) 
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d A 
dt 



* (A 

[i a 



2£] 

0e 



4* A 



2 h 

0e 



) + 






( A a fi + A e £ 2 + A E ) 4^ (56) 



lit 

dt 



7 U a 



ii 

3E 



4- A 



lh 

0E 



} + ^2 ( A a f l +A e f 2 + A E> ^ 



0 E 



(59) 



B. Control Equation for Optimal Thrust Angle 
For the optimal thrust angle: 



an 

da 



= 0 = 



0 f 

a 0 



4* A 



e 




(60) 



Equations (57), (58), and (60) are identical to equations (35), (36), and 
(41), using dE = — dt. Thus, the equation for tan a is the same 
as equation (50) : 



tan a = 



2a A X 4* A Z M 

a e — _N_ 

2a A Y 4* AW - D 

a e 



(61) 



Also, since t is not explicit in equation (52), the Hamiltonian 
is a constant: 



A f , 4- A f 2 + ^ p 

H = — — ^ = CONSTANT 



(62) 



C. Transversality Condition 

The transversality condition is now 



Therefore 



t - Hdt 4* A .dxl 

1 do 

[■ 



+ dG = 0 



Hdt + A da + X de + A„dE 
a e b 



+ dt^ = 0 



(63a) 



(63b) 
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Applying the boundary conditions from equations (38) and (39), 



dt , da , de , dE , da., de_ = 0 
o o o o f f 

- Hdtj- + dE^ + dt^ = 0 



Hence, since E^ is unspecified, 



and 



A E f 2A E< t f) ' 0 
since t^ is also unspecified, 



- Hdt^ + dt^ = 0 
which yields H = 1 

Combining equations (62), (64) and (66) 

A a £ l + A e £ 2 + A E 



(63c) 

(63d) 

(64) 

(65) 

( 66 ) 






= 1 



(all t) 



and : 



I = A a f l +A e f 2 ~ 13 = ° (t = t f ) 



(67) 

( 68 ) 



Using equations (67) and (68) the Euler-Lagrange system of three 
equations (57), (58) and (59) can be reduced to the following two equations: 

0f, 0f„ 



d A n 

0 __a = M 



dt 



9a 



dA 

e = 

dt 0e 



d& 



-(a 9 . + A 

\ a d a e d a f 

/ Ifi . 

y a ? e A e de J 



The solution of equations 



da 

dt 

de 

dt 

dE 

dt 



P 



(69) 

(70) 

(71) 

(72) 

(73) 



together with equations (68), (69) and (70) is obtained numerically by choosing 
initial values as follows: 

C o = °> a o ’ e o ’ E o = ° ’ A a (0) ’ x e (0) (7A) 

This requires guessing of values for A^(0) and A (0) since there is 

at present no analytic way of choosing them. They must be chosen so that the 



numerical integration yields the required set of final values as follows: 



a - a 



f ’ 



e - e 



f ’ 



( 



t f 



(75) 
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V. THE MINIMUM E PROBLEM 

As mentioned earlier, the minimum E problem sheds some light 

on the minimum t problem. We now will minimize the final E (number of 

turns going from a,e,E =0 to a. , e„). 

o o o f f 

A. Hamiltonian and Euler-Lagrange Equations 

Using equations (30) and (31) minimize G = E f . The Hamiltonian 
is 



H A a f l + A e f 2 

The Euler-Lagrange equations are 



a _ v u . i . 

dE 0a ^^a0a + ^e0a' 



(76) 



(77) 



d X 



e 



dE 



0 H 
0e 



( A 



0 e 



+ 




(78) 



These are different from the Euler-Lagrange equations for minimum 
t optimization. 

B„ Control Equation for Optimal Thrust Angle 
For the optimum thrust angle: 

JJl . o - A -¥l + x (79) 

0a a^« eda 
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C. Transversality Condition 

Now, writing the transversality equation. 



-HdE + 



*1 f 



X .dX. 




+ dG = 0 



(80a) 



where 



dG = dE f 



(80b) 



and expanding equation (80a) 



o o o o 

-Hd£ - HdE £ + 1 d/ + )i d/ r + X di + 
o f a o af eo 

o 

X e d 4 + d Ef - o 



(80c) 



There fore 



- HdE £ + dE £ = 0 
f f 



(80d) 



which yields 



H 



f 



X a f l + X e f 2 (at E E f ) 



1 



(80e) 
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VI. TRAJECTORY ANALYSIS 

A. Comparison of Minimum t and Minimum E Trajectories 

Equation (80e) is significant; it means that any pair of initial 

A's (A 3q , A e ) will automatically generate minimum E trajectories 

terminating at any a^ , e f ’ anc ^ ^f * This is so because the condition 

of equation (80e) can always be met simply by scaling the A 's . Note that 

equations (77) and (78) are both linear in the A 's and that the optimal 

thust angle control equation (79) relies only on the ratio of A /A and 

a e 

not on their separate magnitudes. Automatic satisfaction of equation (80e) 

makes it easy to calculate minimum E trajectories and to construct a 

graphical representation of minimum E transfers on an e/a plane. 

Consequently, many transfer trajectories are generated which optimize 

E , and also some which optimize t . These optimal orbital transfers can 

be seen in Figure 3, which is a graph of e vs a/a^ . It shows both minimum t 

and minimum E trajectories which are generated from the same initial orbit 

around Venus. These particular transfers concern a problem incorporating 

solar-electric orbiters around the planet Venus. (The thrust acceleration 

-3 2 

used in the calculations is 5x10 m/s ). The initial gravitational acceler- 

2 

ation is about 0.8 m/s ; thus the basic assumption that A <( <^ g is satisfied. 

In an actual case the value of thrust acceleration would be about 

-3 2 -3 

1/25 of 5x10 m/s . However, the value of 5x10 is used for numerical 

calculations in order to speed up the trajectory and computer times, and thus 

reduce computer costs. Thus all calculated times must be scaled up by a 

factor of 25 to achieve a more realistic mission time. 
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Figure 4 shows minimum E trajectories only, with lines of 
constant t and E superimposed. Examination of the minimum t and 
minimum E trajectories readily discloses that orbital transfers from one 
point on the graph to another are accomplished in the least time by following 
the minimum t path and in the least total revolutions by following the 
minimum E path. Figure 5 shows the minimum t trajectories superimposed 
on the constant t and E lines of Figure 4. The minimum E trajectories 
are omitted for clarity’s sake. 

The actual calculated orbital elements oscillate over each turn, 

as illustrated by curve 1A in Figure 5. For simplicity, the remaining curves 

are constructed by drawing a smooth curve through the values of the orbital 

elements of each turn at E = 2n?r (n = 0, 1, 2, ...). Curve 1A shows the exact 

path of a minimum t trajectory. It is obviously an oscillatory type motion, 

which is much more pronounced here than it would normally be due to the 

increased thrust acceleration (25 times greater than normal). With the actual 

-4 2 

thrust of about 2x10 m/s the oscillations would adhere much more closely 
to the average path depicted by the heavy dotted lines - (curve IB) in 
Figure 5. The minimum t trajectories, which are also minimum fuel trajectories 
here, are very similar to the minimum fuel trajectories calculated with variable 
thrust by Edelbaum. 

B. Rotation of Line of Apsides 

When considering optimal orbital transfers between coaxial orbits, one 
wishes to maintain the same line of apsides. Since this program is designed 

See Section VIII. 
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to change only a and e , no change in the line of apsides is expected. 

A close watch is kept of the change in apsidal line angle throughout all 
orbital transfers and it is verified through calculated results that such 
changes are small and relatively insignificant. In some instances a change of 
up to one degree in the line of apsides is detected over a transfer involving 
up to about 50 revolutions, i.e. , one degree change in 18,000 degrees of 
revolution. 

Several points along curves IB and 2 in Figure 5 indicated along 
with the flight time and eccentric anomaly elapsed at each position. 

Discussion of these points follows. 

Point J . Following a minimum time trajectory, the vehicle arrives 
at point J after 3.03 days* time and 37.69 radians of revolution. This 
point is clearly well past the three-day minimum E line and by inter- 
polation it is easy to see that a minimum E trajectory would take about 
3.2 days to transfer to the same orbit. Now, while the minimum time 
trajectory has gotten to the same point in less time than the minimum E 
trajectory, it has gone through more degrees of revolution to accomplish 
this. Thus the minimum time trajectory has gone through 37.69 radians of 
revolution; while, by interpolation, it can be seen that the minimum E 
trajectory will reach the same point with only about 36.8 radians of 
revolution. In summary, a minimum time transfer takes 3.03 days in 37.69 
radians, while the minimum E transfer takes about 3.2 days in 36.8 radians. 

Point L. Here the minimum t trajectory has reached che point L 
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in 4.57 days while it takes the minimum E trajectory about 4.8 days. How- 
ever, it takes the minimum t trajectory 50.26 radians of revolution as 
opposed to only about 47 radians of revolution for the minimum E trajectory. 

Point P . This is a very good point for comparison. The minimum t 
flight time is 4.9 days. This is clearly less than the minimum E time, 
which is well past the five day line. The eccentric anomaly elapsed is 
46.9 radians which clearly exceeds the value for the minimum E trajectory 
since point P is distinctly below the minimum E contour line where 



E - 45.23 radians. 



0.9 



o 

o 




FTfnmK 



MINIMUM t AND MINIMUM E TRAJECTORIES 
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TRAJECTORY 
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VII. DISCUSSION OF OSCILLATORY MOTION 

A. Effect of Thrust Acceleration and Lagrangian Multipliers on Oscillatory Motion 

Through numerous computer runs, it was verified that different 
thrust accelerations produce essentially the same mean transfer trajectories, 
providing the Lagrangian multipliers are scaled according to the ratio of 
the accelerations (observe Figure 6). A very low thrust system produces a 
trajectory with minimal oscillatory motion (curve A^) . As the thrust 
acceleration is increased, and the Lagrangian multipliers scaled accordingly, 
the oscillatory motion of the transfer trajectory becomes much more 
pronounced (curve A^) , while the mean trajectory remains the same as with the 
lower thrust acceleration. 

It has also been verified (Figure 7) that minute changes in the initial 
Lagrangian multipliers result in uniquely different transfer trajectories 
which do not overlap, as long as the same thrust acceleration is used. The 
transfer trajectories can overlap if both the Lagrangian multipliers and the 
thrust acceleration are changed (Figure 8). 

Figure 9 shows orthogonality between minimum E trajectories and 
constant E path which must exist for optimal solutions. This is a good 
check on the problem solution. In this case, each point was plotted to show 
^the detailed angular orientation of the two curves at the point of inter- 
section. The smoothed-in path of trajectory P in Figure 4 does not appear 
to show orthogonality, but close observation, i.e., blow-up in Figure 9, proves 
it to be so. 

B. Appearance of Negative Eccentricities 



An interesting and perplexing problem was encountered while examining 
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the oscillatory points along trajectories at very small eccentricities. 

The problem is one which involves explaining the presence of a few negative 
eccentricity values which appear during the calculation of transfers of low 
eccentricity. Figures 10 and 11 illustrate this point. In the blow-up of 
area I in Figure 11, points A, B, C and D all lie below the axis and are 
thus negative eccentricities (of very low magnitude) and are printed out as 
such by the computer. 

These negative points appear to be in line with the normal cycle 
of the oscillatory motion. Another point to note is that none of the 
smoothed trajectories ever go below zero eccentricity; that is to say, at 
the end of every cycle (i.e., for every 2mr of eccentric anomaly) e 
is always greater than or equal to zero. However, to date there is no clear 
meaning accorded to these negative eccentricities. Further analysis is 
necessary to determine the real significance of this factor. 

C. Oscillatory Motion and the £ Function 

Also connected with the oscillatory motion which accompanies all 
spiral trajectories is the £ function as defined in equation (40). Through 
many minimum t program runs, it was found that the £ function never stays 
at zero but always oscillates through zero once every revolution, usually at 
every 2nn of eccentric anomaly. The smoothed £ function is zero, but the 
minimum t formulation states that the £ function must equal zero at the 
terminal point in order to represent an optimal transfer. Therefore, the 
transfers are exactly optimal once every revolution when the £ function is 
zero and are near optimal in between. The minimum E transfers are exactly 
optimal at all points along the trajectory since there is no £ function 
requirement in that formulation. 
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A, > A 



COMPARISION OF TRAJECTORY OSCILLATORY MOTION 
FOR DIFFERENT VALUES OF THRUST ACCELERATION 



FIGURE 6 
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COMPARISON OF TRAJECTORY OSCILLATORY MOTION 
FOR DIFFERENT VALUES OF LAGRANGIAN MULTIPLIERS 



FIGURE 7 
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COMPARISON OF TRAJECTORIES FOR DIFFERENT 
VALUES OF THRUST ACCELERATION AND LAGRANG- 
I AN MULTIPLIERS 



FIGURE 8 
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BLOW-UP OF MINIMUM E TRAJECTORY P SHOWING 
ORTHOGONALITY BETWEEN THE TRAJECTORY AND THE 
CONSTANT E PATH WHICH MUST EXIST FOR OPTIMAL SOLUTIONS, 

FIGURE 9 
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OCCURENCE OF NEGATIVE ECCENTRICITIES 



FICUKt: 10 
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VIII. COMPUTATIONAL COMPARISON WITH EDELBAUM'S VARIABLE THRUST CALCULATIONS 
A. Examining Thrust Angle Schemes for Different Trajectories 



limit of best possible performance with an optimal fixed-thrust system. 

This is so because of the extra degree of freedom added to the control system 
when the variable thrust is incorporated, making it a more efficient system 
than the fixed-thrust one. Thus, any optimal performance by the fixed-thrust 
system can, at best, be only equal to the optimal variable thrust system 
and, more often than not, it will be less than the variable thrust system’s 
performance . 



for optimal variable thrust trajectories using an averaging technique and are 
used here for comparison with the fixed-thrust performance: 



An optimal variable thrust system can be utilized to measure the 



The foil wing equations were derived by Edelbaum (References 1, 2, 3) 




(81) 




(82) 



2 

A = the averaged squared variable thrust acceleration over a 
turn = a constant 



(83) 




t 



f 



o 



(84) 
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P = POWER = j | m j v^. 2 



Equation (81) applies to large changes in (a, e) ; i.e., more 

than one revolution. Equation (82) is a special case of (81) and applies 

only to changes within one turn. 

To compare the two propulsion systems, V and aV for fixed 

c c 

thrust are both defined as: 

V = AV = A x t (At) 
c c 

As can be seen in Table 1, the fixed thrust V and aV are 

c c 

always larger than the corresponding variable thrust V or A V where: 

c c 

c = (AV) fixed 
1*" (AV) variable 

Note that is always greater than 1. This indicates a higher cost for 

the fixed thrust system as expected. 

As might be anticipated in cases involving near-circular orbits 
the thrust magnitude of the variable thruster is almost constant, yielding 
the same case as the fixed thrust system. Here the value for is very 

close to one. When transfers are conducted with more elliptical orbits, 
the variable thrust system thrust magnitude scheme no longer resembles 
the fixed thrust case, and the cost is then less, making C ^ take on larger 



(85a) 

(85b) 

( 86 ) 

(87) 



values . 
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B. Maximum A a and Ae Steering programs 

It is interesting and instructive to see how the thrust angle varies 
throughout each transfer. Though only the angle varies (and not the magnitude 
of the thrust) in fixed thrust transfers the angle program agrees closely to 
the one calculated by Edelbaum (See Figure 12). The thrust at points A through 
E coincide in both direction and magnitude. 

Examining single turn characteristics, if e is changing more 
rapidly than a , then the second term on the right-hand side of equation (82) 
will dominate the first term. Similarly, if a is changing more rapidly than 
e , then the first term dominates the second term. Figure 12 shows the 
thrust schemes for maximum A a/a and maximum Ae programs. The last two 
columns in Table 1 show that when the predominant change is in a , then the 
thrust scheme is an "a" thrust scheme and vice versa. 

The fact that Table 1 indicates close but slightly inferior 
performance for the fixed thrust when compared with the variable thrust lends 
considerable support to the validity of both sets of calculations. 
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MAXIMUM Ae 



FIXED 

VARIABLE 



COMPARISON OF FIXED THRUST WITH VARIABLE 
THRUST FOR TRANSFERS WHICH START WITH CIRCULAR ORBITS 

L l> >M: l I .♦ 



Comparison of Fixed-Thrust Acceleration with 
Variable-Thrust Acceleration Performance 
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IX. DETAILED ANALYSIS OF THE THRUST ANGLE SCHEME FOR A MINIMUM E 
TRAJECTORY 

The minimum E trajectory L from Figure 4 was selected for 
examination of its thrust angle scheme. The optimal, computed thrust angle 
was recorded and plotted versus eccentric anomaly (see Figure 13) . The 
trust angle positions were recorded every 90 degrees of eccentric anomaly 
for every third revolution. 

Figure 4 shows that for the first 2.6 days the thrust program 
changed primarily the semi-major axis. In this time period a doubled from 
2 to 4x10^ meters while the eccentricity increased very slowly from 0 to 
about 0.03. This part of the transfer is depicted in Figure 14 where the 
original orbit is 0^ , and after 2.6 days the trajectory has reached 
orbit 0^ . Thus the program starts out by increasing the semi-major axis 
with little change in eccentricity; it uses essentially a maximum Aa thrust 
angle scheme with mostly tangential thrust. This can be seen by observing 
the thrust angle scheme for revolutions one and four in Figure 13. 

By observing trajectory L in Figure 4 one can see that between 
t = 2.6 and t = about 4.6 (days) the thrust program changes to a maximum 
Ae trajectory, while remaining fairly static in its semi-major axis. 

Figure 15 depicts this leg of the transfer from orbit 0^ to orbit 0^ . 

These orbits were chosen at appropriate points to illustrate the computed 
trajectory. A drastic change in the thrust angle scheme between revolutions 
four and seven can be seen in Figure 13. 

Figure 16 shows the transfer between t = 4.6 and t “ about 11 days. 
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The semi-major axis a decreases steadily while e increases. The engine 
is cut off when the trajectory reaches orbit 0^ . At this point e = 0.998 
and a = 1.24x10^ meters. It is interesting to note that between t = 4.6 
and 11 days, as the program begins to change a more rapidly, the thrust 
angle pattern (between revolutions 10 and 27) begins to revert back again 
toward a tangential type thrust scheme (Figure 13). 



REVS E (radians) 
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THRUST ANGLE SCHEME FOR MINIMUM E TRAJECTORY L 

FIGURE 13 




STARTING FROM ORBIT 0, THE SPIRALING PROCESS HAS 
REACHED ORBIT 0 ATt = 2.6 DAYS. 



FIRST LEG OF MINIMUM E TRAJECTORY L 
SHOWING MAXIMUM Aa STEERING SCHEME 



FIGURE 14 
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AT t ■ 4.6 DAYS THE ORBIT SPIRALING PROCESS HAS 
CARRIED THE VEHICLE TO 0 3 . 



SECOND LEG OF MINIMUM E TRAYECTORY L 
SHOWING MAXIMUM Ae STEERING SCHEME 



FIGURE 15 
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AT 1=11 DAYS, THE VEHICLE HAS BEEN 
TRANSFERRED TO ORBIT 0 4 AT THIS POINT THE ENGINE IS 
CUT OFF AND THE VEHICLE REMAINS IN THE SPECIFIED ORBIT 0 4 • 



FINAL STAGE OF MINIMUM E TRAJECTORY L 
SHOWING ATTAINMENT OF DESIRED FINAL ORBIT (0 4 ) 



FIGURE 16 
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X. MARS ORBITAL TRANSFER 
A. Description of Problem 

In order to gain a practical comparison with single impulse 
chemical retro transfers, realistic values were assigned to solar electric and 
chemical rocket engines to be used in an orbital transfer about the planet 
Mars. An initial elliptical orbit was chosen with large eccentricity and 
large semi-major axis (when compared to the radius of Mars) . The problem 
involved transferring from this initial orbit into a much smaller circular 
orbit with radius of 1.1 that of the radius of Mars. The periapse of the 
initial elliptical orbit is equal to the radius of the desired final circular 
orbit . 

The low-thrust spiral transfer is accomplished by turning on the 
solar electric engine at the initial ellipse periapse, causing the space 
vehicle to slowly spiral into the desired final circular orbit using the 
optimized thrust angle program for fixed thrust. The retro transfer is 
accomplished by using a single chemical rocket retro at the ellipse periapse 
causing the space vehicle to transfer directly into the circular orbit. The 
amount of fuel required for each transfer is then compared. The purpose of the 
spiral transfer is to deliver a larger payload to the final orbit. One must 
then decide if the increase in payload delivered is worth the extra time required 
to make the transfer. Figure 17 shows the two methods of transfer being 
compared, the upper drawing being the spiral transfer and the lower one depicts 
the single chemical retro impulse transfer at the ellipse periapse. 
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INTERMEDIATE 

ORBITS 




SPIRAL TRANSFER 

o 0 = 9. 734764 X lO 7 METERS 

e Q =0.95773909 




CHEMICAL ROCKET RETRO 

mars orbital transfer 

FIGURE 17 



B. Initial Data Used 
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The values used in this comparison were as follows (Reference 8) : 
Solar Electric Engine 

F = 0.59756046 - 

v . = 39 ,464 . 9 m/s 
J 

Chemical Rocket Engine 

v . = 2890.49 m/s 
J 



Both Engines 



M - - = initial mass in ellipse - 2049.5093 kg 
ell 



radius of Mars: 3400 km 



,13 3 , 2 



k = 4.293x10 m /s' 
Mr ' 



lars 



a - - - 9 . 734764x10 m 
ell 



e = 1 - 



ell 



= 0.95773909 



C. Chemical Rocket Transfer 

The Mathematical Development follows: 

v e = velocity at ellipse periapse 



■ 2k _ k 
r P a 



( 88 ) 



= 4.6389x10 m/s 



v 

c 



AV 



c 



= velocity at final circular orbit 

= \ — = 3,390 m/s 

* a 

=* v - v = 1 , 248 . 9 m/s 
ep c 



(89) 



-48- 



m„ 



mass in final circular orbit after 
chemical retro 



m Cl /m 



1 ell 



e 



-Av/v. = -1248,9/2890.49 

j e 



-.431 

e = .652 



So m = (.652)(2049.5093) = 1,340 kgs 



Am c = 2049.5093 
- 1340.0000 
709.5093 kgs 

Thus, the direct chemical rocket insertion requires 709.5093 kg 
of fuel with an initial total vehicle mass of 2049.5093 kg. This means 
that a total mass of 1,340 kg is delivered into final orbit by this method. 
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D. Optimized Fixed-Thrust Acceleration Spiral Transfer 



Next, the same transfer is accomplished using the solar electric 
spiral technique. Assuming the same initial ellipse and total vehicle 
mass, the initial thrust acceleration is determined: 

F -4 2 

A 2 . 92x10 m/s 

m ell 

Although assumed a constant, A will actually increase as the 

-4 2 

vehicle mass decreases due to loss of fuel. A mean value of about 3x10 m/s 
was used in the programs, which implies a total anticipated change of about 
6 percent in thrust acceleration. The time to spiral optimally to e^ , a^ 
was computed to be 128.685 days. To compute the fuel expended while spiraling 
into the final circular orbit (Reference 9): 

F = | m [ v (91) 

l.l _F_ = ,59756046 kg. m/s 2 

I m I v. 39,464.9 m/s 

J 

m = 1.515xl0' 5 kg/s 
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To find the fuel expended while spiraling, Am s , multiple |m 
by the total time t elapsed during the spiraling process: 



Now, 



Am ~ 
s 




t 



t = 128.685 = 1 . 11716xl0 7 sec 
Am s = (1.515xl0‘ 5 kg/s) (1.11716xl0 7 s) 



(92) 



=169 kg* 

E. Results and Comparison of Data 

Thus, the solar electric spiraling process has accomplished the 

same transfer while using only 169 kg of fuel. This yields 



Initial mass 2049.5093 

Fuel - 169.0000 

m 1880.5093 kg 

c 2 

1880.5093 kg of mass is delivered into final circular orbit by the 
solar electric method. Since 709.5093 kg of fuel were required by the chemical 
retro transfer, and only 169 kg were expended during the solar electric spiral, 
the latter method yields a savings of 540.5093 kg or about 1150 pounds. 

These results mean that the desired final circular orbit can be 
obtained almost instantaneously by chemical insertion, but with a cost of 
709.5093 kg of fuel; or, the same orbital transfer can be accomplished over a 
period of 128.685 days with a cost of only 169 kg of fuel. Thus, if one chose 
the slower spiral technique, the amount of payload delivered to final orbit 
could be increased by about 1150 pounds. The advantages and disadvantages of 

* 

Note that for constant thrust magnitude, the total increase in thrust 
acceleration due to fuel depletion = 9 percent. 
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each method must be weighed along with the mission objectives to arrive at 
a decision over which method should be used. It would seem that the spiral 
technique as described would be far better suited to an unmanned observational 
missions than the chemical method (Reference 10) . It would seem impractical 
to use the slow spiral technique on any manned flights, regardless of the fuel 
saved, but this does not preclude the use of solar electric engines enroute 
to various destinations as a thruster boost. 
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XI. CONCLUSIONS AND RECOMMENDATIONS 

The numerical solution of the optimum fixed-thrust acceleration 
elliptic spiral technique has been used to generate many different planeto- 
centric transfer trajectories. Figure 2 represents an infinite number of 
minimum time and minimum eccentric anomaly transfers. It should be noted that 
any transfer generated can also be reversed by appropriate manipulation of 
the Lagrangian multipliers. Optimized transfers can be conducted from any 
point along a generated trajectory to another point on that trajectory. This 
is done by starting the program with the Lagrangian multipliers already 
printed out at the desired point along the trajectory. 

At present there is no analytic method for choosing the Lagrangian 
multipliers for particular transfers. By generating many transfer trajectories, 
one can accumulate many cases, but the initial process involved in choosing 
the Lagrangian multipliers remains guesswork and experimentation. Possibly 
by examining more closely the transfers produced by the numerical solutions 
used herein one might draw some clues for the analytic solution. The use of 
reiteration techniques to converge initial guesses to final correct values 
was not explored in this thesis, but is a well-established art (Reference 11). 

As mentioned in an earlier section, all transfers were designed to 
be coplanar, and co-axial. As expected, virtually no rotation of the line of 
apsides was observed. A problem noted and recommended for further study was 
the occasional production of very small negative eccentricities along certain 
transfer trajectories. It only occurs in transfers involving circular or 



near-circular orbits and then only in transfers displaying maximum thrust 
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angle schemes. 

The hypothetical Mars orbital transfer conducted under realistic 
conditions yields some interesting facts for consideration. The optimized 
spiral technique proves to be able to deliver up to 1150 pounds more into the 
desired final orbit than the single impulse chemical retro method, based upon 
a decrease in the amount of fuel required to accomplish the transfer. However, 
the spiral transfer requires 128 days to make the transfer, while the chemical 
method works almost instantaneously. 

Another possibility would be to assume that no more payload is 
required on this particular flight, in which case the extra 1150 pounds 
could be trimmed from the initial mission weight. This would greatly reduce 
the fuel required to launch the vehicle from earth orbit; or, one could use 
the same amount of booster fuel in which case the vehicle would reach the 
target planet much faster than before and would thereby decrease the overall 
flight time. One must then look at all of the factors involved in particular 
missions to determine which method or combination of methods might be most 
useful . 

While this thesis has addressed itself only to the use of optimized 
fixed-thrust orbital transfers, it is quite conceivable, and possibly even 
more advantageous to consider using the electrostatic engine in heliocentric 
space enroute to the target planet as well as after planetocentr ic space has 
been entered. Any judgment on how to interface a heliocentric interplanetary 
trajectory with planetary orbit insertion will have to be based upon further 
study which must involve optimizing the entire flight from earth orbit to the 
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desired final orbit about the target body. Figure 18 illustrates this type 
of problem. This involves the incoming parabolic trajectory, chemical retro- 
insertion into a very large eccentric ellipse, followed by the optimized 
spiraling transfer technique until the desired final orbit is obtained; in 
this case the final orbit is a small circular one. 

It seems obvious that the slow spiral method of transferral would 
not be desirable for manned flights because of their long duration, but 
there are some instances where such an alternative might be preferred. For 
unmanned exploratory flights, where time may not be the most significant 
factor, this method offers an excellent chance to study the upper atmosphere 
of a given planet and to make planetary observations in a low, systematic 
manner. Starting in a large elliptical orbit, a spiral trajectory which 
steadily decreases in altitude until a low circular orbit is achieved might 
be desirable to record radiation, magnetic, and electrical fields at 
systematically decreasing altitudes, along with other upper level atmospheric 
conditions, planetary observations, etc. Orbital changes conducted by large 
impulse chemical retros make it impossible to record such long-time data as 
the transfers occur so quickly. There are many factors to consider when 
weighing the merits and drawbacks of optimized electrostatic propulsion 
(Reference 10) . 

The surface has only been scratched in this field. Improvements 
lie in many directions. For example, Handelsman (Reference 12) has put forth 
a formulation for the introduction of a switching function into the problem 
which would optimize transfer trajectories for minimum fuel and fixed-time 
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FINAL ORBIT 
(AFTER SPIRALING) 




MARS RENDEZVOUS -SPIRAL TRANSFER 



FIGURE 18 
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restraints. This would allow part of the trajectory to consist of coasting 
arcs in which the engine would be shut off by the switching function in order 
to comply with the fixed-time restraint. 

Based upon the work done to date the following recommendations for 
further study are made. 

1. Investigation of the occasional occurrence of negative 
eccentricities in certain transfer trajectories. 

2. Use of a thrust acceleration which, for fixed-thrust magnitude, 
increases as the vehicle mass decreases due to fuel expenditure. 

3. Optimization of the entire flight trajectory including inter- 
facing with the rendezvous and final orbit insertion in order to find the best 
place to use electrostatic thrusting in the flight plan. 

4. Study the solutions generated with a switching function for 
fixed-time problems and compare with solutions of the type covered in this 
thes is . 

5. Addition of other orbital changes, such as rotation of apsidal 
line, non-planar transfers, etc. 

6. Investigation of an analytic solution for a method of selecting 
the correct initial adjoint variables required for particular transfers 
desired. This would eliminate the guesswork now required at the outset of 
every problem. If this is not successful, iteration techniques may be 



necessary. 
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APPENDIX A 

Description of Computer Program 

This is a description of the computer program used for the numerical 
solution of optimal (minimum time and minimum E) planetocentric spiral 
trajectories using fixed power, fixed thrust propulsion, where the thrust 
acceleration is much less than the local gravitational acceleration. The 
program is written in Fortran IV computer language and incorporates a 
straightforward Runge-Kutta type integration process. It was designed for 
use on the 360/30 computer system, but can easily be adapted for use on 
other systems. Double precision was used in all program calculations 
(Reference 13). 

In this appendix the important program input and output parameters 
are described, followed by the program itself and the program flow chart. 

Program Input Parameters 

In order to produce a trajectory for optimal orbital transfer, 

certain initial quantities must be fed into the computer. These can be 

divided basically into three sets of values; the initial orbital conditions, 

the two Lagrangian variables (A >A )> ar *d the thruster conditions. 

a e 

The initial orbital conditions include: 

a (initial semi-major axis) 
o 

e (initial eccentricity) 
o 

k (mass constant) 

t =0 (transfer time) 
o 

E =0 (eccentric anomaly) 
o 

w = 0 (angular position of line of apsides) 
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The thruster conditions include: 

A (thrust acceleration) 

Ml (number of revolutions desired) 
M2 (number of steps per revolution) 
The initial Lagrangian variables are: 

A (°)s 

a a~ 



^(0)5A o 

e e. 



The accuracy is determined by M2 (the number of steps per revolution 
desired). In other words, if M2 = 50, then the program will compute and 
print out the spacecraft’s coordinates 50 times during each revolution, giving 
the exact description of the elliptical orbit which would be obtained were 
the engine to be cut off at any one of those points. However, we have found 
through experimentation that a value of M2 = 10 is sufficiently accurate and 
generates trajectories more quickly. 

For every initial set of Lagrangian variables, a different transfer 
trajectory will be generated. The type of transfers generated are indicated 
in Figure 4. 



Computer Output Parameters 

At each point in the transfer trajectory a number of instantaneous 
quantities are printed out. These include the total step number, which, if 
divided by the number of steps per revolution, yields the number of revolutions 
to that point (a total step number of 156 would mean 15.6 revolutions at 
ten steps per revolution); the eccentric anomaly in radians; the thrust angle 
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in degrees; the £ function, which must go through zero for an optimal transfer 

to take place; the eccentricity; the change in eccentricity per radian of 

eccentric anomaly; the semi-major axis; the change in semi-major axis per 

radian of eccentric anomaly; the Lagrangian variable A and its change per 

a 

radian; the Lagrangian variable A ^ and its change per radian; the time 
elapsed since thrust initiation; the number of seconds of time per radian of 
eccentric anaomaly; angle of the line of apsides; the change in the line of 
apsides per radian of revolution. 




. 
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Cdola) MAIN 



Program Listing 



DOUBLE PRECISION W , X , Y , Z , A ,GMU , N , D , H ,SIGN 
COMMON W, X, Y,Z, A,GMU, N, D, H, SIGN 

DOUBLE PRECISION HI (4) , H2 (4) , H3 (4) , Q (7) , QQ (7) , DQ (7) , DELQ (7) ,D1,D2 
1 T, TT, PR (7) 

NAMELIST /INPUT/ M 1 , M 2, A, GMU, Q , S IGN 
DATA Q / 7*0. DO / 

SIGN = 1 . DO 
DO 10 M = 1, 4 
H 2 ( M ) = 1 . DO 
HI ( M ) = 0.5D0 
10 H 3 ( M ) = 0.5 DO 
HI (3) = 1 . DO 

113 (4) = 1 . DO 

H3 (1) = 0. DO 
H2 (1) = H 2 ( 1 ) / 6. DO 

H2 (4) '= H2 (1) 

112 (2) = H2 (2) / 3. DO 
H2 (3) = H2 (2) 

20 READ (1 , INPUT, END=200) 

IP( Ml ) 25, 200, 25 

25 WRITE (3, 40) Ml, M2, A, G MU , SIGN 

40 FORMAT ( 1 H 1 , 15, ' REVOLUTIONS', 15,' STEPS PER REV. A = ' , 

1 1PD15.8, ' MU = ' , D 1 5 . 8 ,' SIGH =',D15.8 ) 

WRITE (3, INPUT) 

T = 0 . DO 

K M AX = IA BS ( Ml * M2 ) 

0(7) = 0(1) 

X = M2 

DT = 6.2831853071795865 / X 
KOUNT = 0 

100 TT = T 

DO 101 M = 1, 6 
PR ( M + 1 ) = Q(M) 

QQ (M) = Q (H) 

101 DELQ (M) = 0 . DO 
DELQ (7) = 0. DO 
QQ (7) = Q (7) 

PR ( 1 ) = T * 57.2957795 
PR (6) = PR (6) / 86400. DO 
PR (7) = 0 . DO 

IF (Q (6) .NE. 0. DO.OR.Q (7) . NE. 0. DO) PR (7) = DAT A N2 (Q (6 ) , Q (7 ) ) 

PR (7) = PR (7) * 57.2957795 

W RITE (3 , 50) KOUNT, T, PR 

50 FORMAT ( 1H0, I6,1P8D14.6 ) 

CALL DERI V ( T, Q, DQ ) 

PR ( 1 ) = DAT AN 2 ( N , D) * 57.2957795 
PR (2) = DQ ( 1 ) *Q ( 3) + DQ (2) *Q ( 4) - DQ (5) 

WRITE (3, 51) PR ( 1) , PR (2) , (DQ (M) , M= 1 , 6) 

51 FORMAT ( 1H , OX, 1P8D14.6 ) 

KOUNT = KOUNT + 1 

K = 1 

GO TO 103 

102 K = 1 + K 

T = TT + DT * H3(K) 
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CALL DE RI V ( T, Q, DQ ) 

103 D 1 = ill (K) * D7 

D2 = 112 (K) * DT 

DO 104 M = 1, 7 

DELO (M) = DELQ(K) ♦ D2 * DQ (M) 

104 Q(M) = QQ(M) + DO (M) * D 1 

I F ( K - 4 ) 102, 105, 102 

105 DO 1 0r> M = 1, 7 

106 0(E) = QO («) ♦ DELQ (M ) 

IF ( KOUNT - K MAX ) 100, 1 07, 100 

107 CALL DERIV( T, 0, DQ ) 

DO 120 8=1,6 

120 PR (M + 1) =Q (M) 

PR ( 1 ) = T * 57.2957795 
P li ( 6 ) •= PR (6) / 86400. DO 

P R ( 7 ) = 0 . DO 

IF (0 (6) .N E.O. DO.OR.Q (7) . NE.O. DO) PH (7) = D ATAN 2 (Q (6 ) , Q (7 ) ) 

PR (7) = PH ( 7 ) * 57.2957795 

W RITE (3,5 0) KOUNT, T, PR 

PR ( 1 ) - DATA N 2 ( N , D) ♦ 57.2957795 

PR(2) = DQ ( 1 ) *Q(3) + DQ (2 ) *0 (4) - DQ (5) 

WRITE(3,51) PR ( 1 ) , PR (2) , (DQ (M) , M=1 ,6) 

GO TO 20 
200 RETURN 
END 
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DEE IV 

IMPLICIT PEAL*8 ( A-H, 0-Z ) 

SUBROUTINE DERIV ( E, Q, DQ ) 

DOUBLE PRECISION W , X, Y , Z, A, GN U , N , D , H , S IGN 
COMMON W , X, Y, Z, A ,GMU ,N, D, H, SIGN 

DOUBLE PRECISION E , Q ( 4) , DQ (4) , DSI N , DCOS , DSQRT, DABS, SINE, COSE, 
1 COS2E, T1 , T2, T 3 
SQHTAU = DSQRT ( Q(2) / GMU ) 

SINE = DS IN ( E ) 

COSE = DCOS ( F. ) 

COS2E = COSE * COSE 
Z = 1.D0 - Q (1) * Q (1) 

Y = DSQRT ( D A BS ( Z ) ) 

Z = Z * SINE 
X = 0(1) * SINE 

W = Y * (2. DO * COS F. - Q ( 1 ) * (C0S2E+ 1.D0)) 

N = 2. DO * Q(2) * Q (4 ) * X + Q ( 3 ) * Z 

D = 2. DO * Q(2) * Q (4) * Y «■ Q (3) * W 

I! = DSQRT ( N*N + D *D ) * SIGN 

T 1 = X *N + Y * D 
T2 = Z*N + W * D 
T 3 = A * Q(2) / GMU / II 
DQ (2) = 2. DO * T3 * 0(2) * Q (2) * T1 

DQ (1) = T3 * Q(2) * T2 

DQ (5) = SQRTAU * (1.D0 - Q(1) * COSE) 

DQ ( 4 ) = -2.D0 * T3 * (3. DO * Q (2) * Q (4 ) * T1 ♦ Q (3) * T2 ) 

1 + 1.5D0 * DQ (5) 

DQ (5) = DQ (5) * Q (2) 

DQ (3 ) = -T3*Q (2) * (2. D0*Q(2) *Q (4) * (N*SINE - Q(1)*D/Y) - 

1 Q ( 3) * (2. DO* X*N + (Q (1) *W/ (Y*Y) + Y*(1.D0 + C0S2E))*D)) 

2 - SQRTAU * Q (2) * COSE 

DQ (6) =T3*Q(2) *(Z*D/(Y*Y)* (2.D0-Q (1) * (COSE + Q (1) ) ) - 

1 Y*N * (COSE - Q (1) ) ) 

T 3 = 0. DO 

IF (Q (6) .NE.O. DO.OR.Q (7) .NE.0.D0)T3 = DATAN2 (Q (6) , Q (7) ) 

T 1 = DSIN (T 3) 

T2 = DCOS (T3) 

DQ (7) = T2*DQ(1) - T1 * DQ (6) 

DQ (6) = T2*DQ (6) - T1*DQ(1) 

RETURN 

END 
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Program Flow Chart 
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SUBROUTINE DF.RIV (E,Q,DQ> 
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APPENDIX D 

Typical Computer Readout 

The following page represents a computer read-out of a particular 
transfer trajectory generated by the minimum time program previously 
described. It was generated at a step size of ten per revolution, i.e., 
step ten is the end of the first revolution. Read across from the step 
desired. The meaning of each column is indicated at its head. 

When reading across one will notice two lines of figures. In 
each case the second figure represents the change in value per radian of 
eccentric anomaly as the trajectory generates, with the exception that the 
second figure in the second column is the £ function which must go through 
zero for optimality. On this page of sample read-out the £ function does not 
yet go through zero. But later on in the trajectory it begins to oscillate 
through zero with every revolution. 
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